PLS-DA of P. acuta gene expression

Note: Open Rproj first, then script. To easily use relative paths, click the down button next to knit and then click “Knit Directory –> Project Directory”. This should make loading and saving files much easier.

1. Load packages and input data

Load packages

library(edgeR, quietly = TRUE) #edgeR-v3.30.3
library(vegan, quietly = TRUE)
library(Dune, quietly = TRUE)
library(ggplot2, quietly = TRUE) #ggplot2-v3.3.5
library(tidyverse, quietly = TRUE) #tidyverse-v1.3.1
library(Rmisc, quietly = TRUE)
library(mixOmics, quietly = TRUE)
#library(ggridges)
#library(hrbrthemes)
#library(viridis)

Load the input file containing the treatment information

#treatment information
treatmentinfo <- read.csv("Sample_Info/samples_Pacuta.annotations.txt", header = TRUE, sep = "\t", fileEncoding="UTF-8-BOM") #read in file
top3groups <- c("Group2", "Group3", "Group6")
treatmentinfo <- filter(treatmentinfo, group%in%top3groups)
table(treatmentinfo$group)

Group2 Group3 Group6 
    34     24     27 

Load the input file containing the gene count matrix

#gene count matrix
gcount <- as.data.frame(read_delim("Genome_Info/Pocillopora_acuta_KBHIv2.gentrome.fa.gz.salmon.numreads.matrix", delim = "\t", col_names = TRUE, show_col_types = FALSE), fileEncoding="UTF-8-BOM") #read in file
rownames(gcount) <- gcount$Name #makes "Name" the rowname
gcount <- gcount[,-c(1)] #drops the "Name" column
gcount <- round(gcount) #round 
dim(gcount); head(gcount)[,1:3] #view dataset attributes

[1] 33259   119

gcount <- gcount[,treatmentinfo$sample]

Determine library size

libSize.df <- data.frame(libSize=colSums(gcount))

Make DGE object

DGEdat <- DGEList(counts=as.matrix(gcount), samples=treatmentinfo,
                  group=treatmentinfo$temp)
dim(DGEdat$counts)

[1] 33259    85

2. Pre-filtering

lib.sizes.cpm <- colSums(gcount)/1000000
lib.sizes.cpm <- tibble(sample=colnames(gcount), lib.size=lib.sizes.cpm)

keep <- rowSums(cpm(gcount) > 3.33) >= 2
table(keep)

keep
FALSE  TRUE 
12211 21048 

DGEdat <- DGEdat[keep, , keep.lib.sizes=FALSE]

3. Data normalization

DGEdat <- calcNormFactors(DGEdat)
DGEdat$samples

4. Plot global gene expression

Log transform the counts matrix for the next plots

DGEdat.cpm <- DGEdat #make a copy the edgeR dataset
DGEdat.cpm$counts <- cpm(DGEdat.cpm$counts, log=TRUE, prior.count=5) #log transform the copy for the next plots

Run a principle coordinates analysis, all samples

d <- data.frame(DGEdat.cpm$samples, name = colnames(DGEdat.cpm$counts), TP_temp=paste(DGEdat.cpm$samples$timename, DGEdat.cpm$samples$temp, sep="_")) #make a dataframe containing all plotting info
#d <- filter(d, temp=="Hot")
d$timename <- factor(d$timename, levels = c("0_hour", "6_hour", "12_hour", "30_hour", "1_week", "2_week", "4_week", "6_week", "8_week", "12_week"))
ncol(DGEdat.cpm$counts[,d$sample])

[1] 85

# Partial Least Squares Analysis
mixOmic_plsda <- plsda(t(DGEdat.cpm$counts[,d$sample]), d$timename, ncomp = 85)
plotVar(mixOmic_plsda, plot = F)

plsda_plot <- biplot(mixOmic_plsda, var.axes = TRUE, cutoff = .8, ind.names = F); plsda_plot

plotIndiv(mixOmic_plsda, ind.names = F, ellipse = TRUE, legend = TRUE)

Run a principle coordinates analysis, Hot samples

d <- data.frame(DGEdat.cpm$samples, name = colnames(DGEdat.cpm$counts), TP_temp=paste(DGEdat.cpm$samples$timename, DGEdat.cpm$samples$temp, sep="_")) #make a dataframe containing all plotting info
d$timename <- factor(d$timename, levels = c("0_hour", "6_hour", "12_hour", "30_hour", "1_week", "2_week", "4_week", "6_week", "8_week", "12_week"))
d <- filter(d, temp=="Hot")
ncol(DGEdat.cpm$counts[,d$sample])

[1] 44

# Partial Least Squares Analysis
mixOmic_plsda <- plsda(t(DGEdat.cpm$counts[,d$sample]), d$TP_temp, ncomp = 44)
plotVar(mixOmic_plsda, plot = F)

plsda_plot <- biplot(mixOmic_plsda, var.axes = TRUE, cutoff = .8, ind.names = F); plsda_plot

plotIndiv(mixOmic_plsda, ind.names = F, ellipse = TRUE, legend = TRUE)

Run a principle coordinates analysis, Amb samples

d <- data.frame(DGEdat.cpm$samples, name = colnames(DGEdat.cpm$counts), TP_temp=paste(DGEdat.cpm$samples$timename, DGEdat.cpm$samples$temp, sep="_")) #make a dataframe containing all plotting info
d$timename <- factor(d$timename, levels = c("0_hour", "6_hour", "12_hour", "30_hour", "1_week", "2_week", "4_week", "6_week", "8_week", "12_week"))
d <- filter(d, temp=="Amb")
ncol(DGEdat.cpm$counts[,d$sample])

[1] 41

# Partial Least Squares Analysis
mixOmic_plsda <- plsda(t(DGEdat.cpm$counts[,d$sample]), d$TP_temp, ncomp = 41)
plotVar(mixOmic_plsda, plot = F)

plsda_plot <- biplot(mixOmic_plsda, var.axes = TRUE, cutoff = .8, ind.names = F); plsda_plot

plotIndiv(mixOmic_plsda, ind.names = F, ellipse = TRUE, legend = TRUE)

Emma script

# set working directory and load necessary packages

library(reshape2)
#library(ggbiplot)
library(broom)  # devtools::install_github("tidymodels/broom")
library(cowplot)
library(ggpubr)
#library(ggfortify)
library(ggrepel)
library(gridExtra)
#library(ggforce)
# set seed
set.seed(54321)

Load datasets

d <- data.frame(DGEdat.cpm$samples, name = colnames(DGEdat.cpm$counts), TP_temp=paste(DGEdat.cpm$samples$timename, DGEdat.cpm$samples$temp, sep="_")) #make a dataframe 
d$timename <- factor(d$timename, levels = c("0_hour", "6_hour", "12_hour", "30_hour", "1_week", "2_week", "4_week", "6_week", "8_week", "12_week"))

#PCA
Mcap.pca.out <- prcomp(t(DGEdat.cpm$counts[,d$sample])) #, center=FALSE, scale=FALSE) #run PCA
M.summary <- summary(Mcap.pca.out); M.summary #view results

Importance of components:
                           PC1     PC2     PC3      PC4      PC5      PC6      PC7      PC8      PC9     PC10     PC11     PC12     PC13     PC14    PC15
Standard deviation     47.8457 39.3804 31.8471 23.44963 19.12449 16.68497 15.27388 14.09774 12.78940 11.69200 11.33499 10.95941 10.53431 10.13147 9.45036
Proportion of Variance  0.2376  0.1609  0.1053  0.05707  0.03796  0.02889  0.02421  0.02063  0.01698  0.01419  0.01333  0.01247  0.01152  0.01065 0.00927
Cumulative Proportion   0.2376  0.3985  0.5038  0.56085  0.59881  0.62770  0.65191  0.67253  0.68951  0.70370  0.71703  0.72950  0.74101  0.75166 0.76093
                          PC16    PC17    PC18    PC19    PC20    PC21    PC22   PC23    PC24    PC25    PC26    PC27    PC28    PC29    PC30    PC31    PC32
Standard deviation     9.07542 8.86547 8.54723 8.36652 8.11764 8.07213 7.72604 7.4751 7.33188 6.98285 6.81872 6.75749 6.59476 6.50150 6.39636 6.34873 6.27945
Proportion of Variance 0.00855 0.00816 0.00758 0.00726 0.00684 0.00676 0.00619 0.0058 0.00558 0.00506 0.00483 0.00474 0.00451 0.00439 0.00425 0.00418 0.00409
Cumulative Proportion  0.76948 0.77764 0.78522 0.79248 0.79932 0.80609 0.81228 0.8181 0.82366 0.82872 0.83354 0.83828 0.84280 0.84718 0.85143 0.85561 0.85970
                          PC33    PC34    PC35    PC36    PC37    PC38    PC39    PC40    PC41    PC42   PC43    PC44    PC45    PC46    PC47    PC48    PC49
Standard deviation     6.15489 6.07861 6.03575 5.97942 5.90674 5.84883 5.79111 5.76026 5.70003 5.66802 5.6388 5.61025 5.55841 5.49090 5.45282 5.42537 5.36089
Proportion of Variance 0.00393 0.00383 0.00378 0.00371 0.00362 0.00355 0.00348 0.00344 0.00337 0.00333 0.0033 0.00327 0.00321 0.00313 0.00309 0.00305 0.00298
Cumulative Proportion  0.86364 0.86747 0.87125 0.87496 0.87858 0.88213 0.88561 0.88906 0.89243 0.89576 0.8991 0.90233 0.90554 0.90866 0.91175 0.91481 0.91779
                          PC50    PC51    PC52    PC53   PC54    PC55    PC56    PC57    PC58    PC59    PC60    PC61    PC62    PC63    PC64   PC65    PC66
Standard deviation     5.33055 5.26422 5.23060 5.20470 5.1901 5.13711 5.12038 5.06910 5.04595 5.01429 4.99416 4.95214 4.93965 4.89052 4.84779 4.8065 4.78845
Proportion of Variance 0.00295 0.00288 0.00284 0.00281 0.0028 0.00274 0.00272 0.00267 0.00264 0.00261 0.00259 0.00255 0.00253 0.00248 0.00244 0.0024 0.00238
Cumulative Proportion  0.92074 0.92361 0.92645 0.92926 0.9321 0.93480 0.93752 0.94019 0.94283 0.94544 0.94803 0.95057 0.95310 0.95559 0.95802 0.9604 0.96280
                          PC67    PC68    PC69    PC70    PC71    PC72    PC73    PC74    PC75    PC76    PC77    PC78    PC79   PC80    PC81    PC82   PC83
Standard deviation     4.76581 4.74457 4.71724 4.67894 4.63921 4.58973 4.55821 4.53962 4.49122 4.47088 4.40430 4.37397 4.34636 4.2770 4.24909 4.21365 4.1665
Proportion of Variance 0.00236 0.00234 0.00231 0.00227 0.00223 0.00219 0.00216 0.00214 0.00209 0.00207 0.00201 0.00199 0.00196 0.0019 0.00187 0.00184 0.0018
Cumulative Proportion  0.96516 0.96749 0.96980 0.97208 0.97431 0.97650 0.97865 0.98079 0.98288 0.98496 0.98697 0.98896 0.99092 0.9928 0.99469 0.99653 0.9983
                          PC84      PC85
Standard deviation     4.00645 5.257e-14
Proportion of Variance 0.00167 0.000e+00
Cumulative Proportion  1.00000 1.000e+00

biplot(Mcap.pca.out) #plot results

Mcapitata.all <- Mcap.pca.out %>%
  augment(d) %>% # add original dataset back in
  group_by(timename, temp) %>%
  mutate(PC5.mean = mean(.fittedPC5),
         PC6.mean = mean(.fittedPC6))

mcap.Hot.segments <- Mcapitata.all %>% subset(temp == "Hot") %>%
  select(timename, temp, PC5.mean, PC6.mean) %>%
  gather(variable, value, -(timename:temp)) %>%
  unite(group, timename, variable) %>% distinct() %>%
  spread(group, value)

mcap.amb.segments <- Mcapitata.all %>% subset(temp == "Amb") %>%
  select(timename, temp, PC5.mean, PC6.mean) %>%
  gather(variable, value, -(timename:temp)) %>%
  unite(group, timename, variable) %>% distinct() %>%
  spread(group, value)